/*!
 * @file randomutil.h
 * @ingroup MathUtil
 *
 * A library for generating random numbers.
 * We currently provide a wrapper around the boost::random generators,
 * but it will be easy to swap this out if need be.  
 * The library is designed to exploit the efficiency of inline
 * functions as much as possible.  Hence, all code is in the header files.
 *
     *  Classes with the names of
     * distributions are mainly for convenience.  They
     * define cdf, inverse cdf, and pdf functions and allow
     * users to bind parameters to an intance of that class.
     * For example, if one wanted a Normal(mu,sigma) distribution,
     * one could specify mu and sigma in the constructor of Normal
     * and subsequent calls to the non-static versions of distribution
     * related functions can be called without the parameters.
     *
     * The exception to this are the classes StdUniform and StdNormal.
     * They represent a single distribution rather than a parameterized
     * family, so there is no need to instantiate them, although you
     * could if you really wanted to.
     *
     * It is assumed that all random number generators can
     * be seeded with any non-zero unsigned long (32 bit adam:check for 64 bit machines).
     * Seed functions and generator constructors
     * should also provide default (void) versions
     * and a version that copies the underlying uniform distribution.
     *
     * Functors and classes not related to the convenience
     * routines should make as few assumptions about
     * the template parameters as possible.  For example,
     * random number generators should only be assumed to
     * support the operator() call.  "Nice" versions of
     * these entities can be wrapped or subclassed in the
     * appropriate distribution class.  They should also
     * provide a uniform factory-type interface to the most common
     * implementations of the generators for that particular distribution.
     *
     *
     *
     *
 *
 */



#ifndef _QT_RANDOMUTIL_H_
#define _QT_RANDOMUTIL_H_

#include <cstdlib>
#include <cmath>
#include <cassert>

#include <boost/random.hpp>


namespace qt {

namespace random {

static const unsigned long default_seed = 11232u;

/*!
 * These typedefs give standard access to candidate base_generator_types.
 *
 */
typedef boost::mt19937               mersenne_twister;
typedef boost::mt11213b             mersenne_twister2;
typedef boost::minstd_rand0         minstd_rand0;
typedef boost::minstd_rand          minstd_rand;
// this one is strange.
//  typedef boost::lagged_fibonacci44497 lagged_fibonacci;
typedef boost::hellekalek1995        inverse_linear;

/*!
   * This class takes a uniform random number generator
   * and a inverse cdf function and produces a generator
   * that gives random numbers with a corresponding distribution.
   */

template<typename InvCdfFuncType,
         typename UniformRandomNumberGenerator>
class InvCdfGenerator  {


private:
    // no default constructor
    InvCdfGenerator();

    // retired constructor...
    //      InvCdfGenerator(InvCdfFuncType invCdf,
    //     unsigned long s=default_seed)
    //: _invCdf(invCdf), _rng(s) {}



public:
    typedef double return_type;

    /*!
     * Constructor needs two arguments.
     *
     *
     * @param
     * @returns
     *
     * History: Still need a valid copy ctor since uniform_01 takes
     * reference to its generator.
     *
     */

    InvCdfGenerator(InvCdfFuncType invCdf,
                    UniformRandomNumberGenerator &rng)
        : _invCdf(invCdf), _rng(rng) {}

    double operator()() {
        return _invCdf(_rng());
    }

private:
    UniformRandomNumberGenerator & _rng;
    InvCdfFuncType _invCdf;

};



/*!
 * This class provides functor typedefs and
 * Holds a base uniform (0,1) random number gererator so that
 * other generators can inherit these methods
 * in a unified way.   This is purely for convenience
 * and uniformity.  This class should not be treated
 * polymorphicaly, i.e., no virtual functions.
 */

template <typename BaseGeneratorType>
class Generator {

private:
    static Generator theGlobalDefaultGenerator;

public:
    /*!
     * This provides reference to a global instance of this random
     * number generator.  Note that there is a separate insatnce
     * associated with
     * each different template specialization.
     */
    static Generator & getDefaultGenerator() {
        return theGlobalDefaultGenerator;
    }

private:
    typedef boost::uniform_01<BaseGeneratorType> uniform_generator_type;

    // no assignment
    Generator & operator=(const Generator &);

    // no compare
    bool operator==(const Generator &);

public:

    typedef double result_type;
    typedef BaseGeneratorType base_generator_type;


    Generator(unsigned long x=default_seed) :
        _baseGen(x), _ugen(new uniform_generator_type(_baseGen)) {}

    Generator(const BaseGeneratorType & gen) :
        _baseGen(gen), _ugen(new uniform_generator_type(_baseGen)) {}


    ~Generator() {
        // this looks kind of strange.  It is done to make
        // copy safe since uniform_01 keeps a reference to
        // a base generator.
        delete _ugen;
        _ugen = NULL;
    }

    void seed(unsigned long x=default_seed) {

        if(0==x) {
            x = default_seed;
        }
        _baseGen.seed(x);
    }

    void seed(BaseGeneratorType bgen) {
        _baseGen.seed(bgen);
    }


    double operator()() {
        return nextUniform();
    }

protected:

    double nextUniform() {
        return (*_ugen)();
    }

private:

    base_generator_type     _baseGen;
    // _gen has a reference to _baseGen, so they must
    // go in and out of scope together.  Hence, the pointer.
    uniform_generator_type * _ugen;

};

// instance of the static member
template <typename BaseGeneratorType>
Generator<BaseGeneratorType>
Generator<BaseGeneratorType>::theGlobalDefaultGenerator(default_seed);


/*! The best all-around generator. */
typedef Generator<mersenne_twister> DefaultGeneratorType;
/*! A slower, but still good generator. */
typedef Generator<inverse_linear> DefaultGenerator2Type;
/*! Fast, but not so good generator */
typedef Generator<minstd_rand> DefaultGenerator3Type;



/*!
 *  This class aggregates a base generator and a distribution
 * class.  This class is again just for convenience and unformity.
 * The call operator can be used to generate random
 * deviates provided the distribution has a invCdf method.
 *
 * This class gains visibility to the generator passed into
 * the constructor.
*/

template <typename GeneratorType, typename DistributionType>
class DistributionGenerator {
public:
    typedef GeneratorType generator_type;
    typedef DistributionType distribution_type;

    typedef double result_type;
    typedef typename GeneratorType::base_generator_type base_generator_type;


    DistributionGenerator
    (generator_type & gen, distribution_type & dist) :
        _generator(gen), _distribution(dist)  {}

    double operator()() {
        return _distribution.invCdf( (_generator)() );
    }

protected:
    generator_type    & _generator;
    distribution_type   _distribution;


};

/*!
 * This class is pairs a Generator object with a WrapperType
 * object that is assumed to hold a reference to a Generator
 * object.  This class owns its Generator.
 */
template <typename GeneratorType, typename WrapperType>
class WrappedGenerator {

public:
    typedef GeneratorType generator_type;

    typedef double result_type;
    typedef typename GeneratorType::base_generator_type base_generator_type;


    explicit WrappedGenerator()
        : _gen(), _wrapper(_gen) {}

    template <typename T>
    explicit WrappedGenerator(T & t)
        : _gen(), _wrapper(_gen,t) {}

    // we may need more constructors here

    void seed(unsigned long x=default_seed) {
        _gen.seed(x);
    }

    result_type operator()() {
        return _wrapper();
    }

private:
    GeneratorType _gen;
    WrapperType _wrapper;
};

}; // namespace qt::random

}; // namespace qt



#endif


